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Topological quantum phases of matter have been a topic of intense interest in contemporary 
condensed matter physics. Extensive efforts are devoted to investigate various exotic properties of 
topological matters including topological insulators, topological superconductors, and topological 
semimetals. For topological insulators, the dissipationless transport via gapless helical edge or 
surface states is supposed to play a defining role, which unfortunately has proved difficult to realize 
in experiments due to inevitable backscattering induced in the sample boundary. Motivated by the 
fundamental connection between topological invariants and the Zak phase, here, we show that the 
non-trivial band topologies of both two and three-dimensional topological insulators, characterized 
by the Chern numbers and the Z 2 invariants, respectively, are directly manifested in the winding 
numbers of the Wannier-Stark ladder (WSL) emerging under an electric field. We use the Floquet 
Green’s function formalism to show that the winding number of the WSL is robust against interband 
interference as well as non-magnetic impurity scattering. 


I. INTRODUCTION 

Formulated in terms of algebraic commutation rela¬ 
tions between operators, quantum mechanics had not 
been usually related with geometry or topology before 
the discovery of the geometric phasd^, or more com¬ 
monly known as the Berry phase. The notion that the 
topological structure of the Berry phase can be used as a 
new “order parameter” distinguishing between different 
quantum phases of matter has triggered an intense out¬ 
burst of research activities in contemporary condensed 
matter physics 3 6 '. 

In two dimensions, where perpendicular spin compo¬ 
nents are conserved, such an order parameter is the spin- 
dependent Chern number, which can be in principle mea¬ 
sured via the spin Hall conductance according to the 
Kubo formula!® G^ in = (Cf — C^)e 2 //i, where C a is the 
Chern number of an occupied band with spin component 
< 7 . Unfortunately, fully spin-filtered measurements are 
very difficult to perform in experiments. An alternative 
is to measure the ballistic two-terminal charge conduc¬ 
tance, which is observed to be quantized approximately 
as G = 2e 2 jh in a two-dimensional (2D) topological in¬ 
sulator (TI) 9 n J agreeing with a theoretical prediction as¬ 
suming that gapless helical edge states generate dissipa¬ 
tionless charge transport 12 . 

In contrast to chiral edge states in the quantum Hall 
effect, however, the helical edge states are inevitably cou¬ 
pled to various backscattering sources induced in the 
sample boundary so that the conductance quantization 
is not exactly protectecP^Sl. This means that the edge 
transport measurement is not an ideal method to reveal 
the bulk band topology of 2D TIs, which remains to be 
protected against non-magnetic impurity scattering as 
well as other perturbations. Given this problem, it is 
beneficial to devise a physical observable directly man¬ 
ifesting the topological order in the bulk without refer¬ 
ence to the non-universal electron dynamics in the sam¬ 


ple boundary. In this context, it is worth mentioning that 
the spin-charge separation in the presence of a 7 r flux de¬ 
fect can be used as a pos sible avenue to reveal the band 
topology in the bulk 2 =d^i 

The situation becomes even more complicated in three 
dimensions, where different spin components are in gen¬ 
eral mixed and thus the spin-dependent Chern number is 
not properly defined. In this situation, proper topologi¬ 
cal order parameters are the Z 2 invariants, (z/o; zff, ^ 2 , ^ 3 ), 
which can fully characterize the three-dimensional (3D) 
band topology if both inversion and time-reversal sym¬ 
metries are present. Unfortunately, no simple transport 
measurement can manifest these 3D topological invari¬ 
ants as directly as the (spin Hall or two-terminal charge) 
conductance in two dimensions. Instead, the strong Z 2 
topological invariant uq (most important for the robust¬ 
ness of a given 3D TI) is predicted to be manifested in 
the quantized magneto-electric effect, where an electric 
field induces a topological contribution to the magneti¬ 
zation 23 . While an actual measurement of the magneto¬ 
electric effect may be too difficult to perform at present, 
the fact that the topological invariant has a direct phys¬ 
ical consequence is important as a matter of principle. 

In practice, the 3D band topology has been inferred 
from the existence of helical surface states, which exhibit 
the spin-momentum locking as a consequence of the non¬ 
trivial topological order 24 25 . However, a problem is that 
the spin-momentum locking by itself is not directly re¬ 
lated with any of the 3D topological invariants. In some 
sense, it is even more beneficial in three dimensions to de¬ 
vise a physical observable directly manifesting the topo¬ 
logical invariants (in addition to the magneto-electric ef¬ 
fect). 

Motivated by the fundamental connection between 
topological invariants and the Zak phase, in this work, 
we show that the non-trivial band topologies of both 2D 
and 3D TIs (characterized by the Chern numbers and the 
Z 2 invariants, respectively) are directly manifested in the 
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energy spectrum of electrons under an electric field via 
the winding number of the Wannier-Stark ladder (WSL). 
The WSL is a set of energy eigenstates of electrons con¬ 
fined in the lattice under an electric field, which are the 
quantized modes of the Bloch oscillation. In contrast 
to a recent interferometric method proposed in optical 
latticed ^ which combines the coherent Bloch oscilla¬ 
tion with the Ramsey interferometry, our spectroscopic 
method can be applied to condensed matter systems, 
where the phase coherence is not guaranteed. Concretely, 
we use the Floquet Green’s function formalism to show 
that the winding number of the WSL is robust against 
interband interference as well as non-magnetic impurity 
scattering. 


Provided that the fully interacting Floquet Green’s 
function is obtained, our method can be applied to any 
strongly correlated systems with the non-trivial band 
topology in order to address various theoretical issues 
such as (i) how far the topological order can persist as a 
function of correlation strength and (ii) what new phases 
of strongly correlated topological matter can emerge at 
sufficiently strong correlation. It is worth mentioning 
that the definitions for topological invariants were pre¬ 
viously extended to general strongly correlated systems 
in equilibrium, i.e., without the electric field 30 32 . Ex¬ 
perimentally, this means that the magneto-electric effect 
(which is directly related with a topological invariant) 
can be used to characterize strongly correlated topo¬ 
logical insulators. In our theory, the non-trivial wind¬ 
ing number of the WSL can play a similar role as the 
magneto-electric effect. 


The rest of the paper is organized as follows. In Sec. [TT| 
we summarize briefly how the 2D and 3D band topolo¬ 
gies are characterized by the Chern numbers and the Z 2 
invariants, respectively. In Sec. |III[ we explain how the 
winding number of the WSL can provide a direct manifes¬ 
tation of topological invariants. Specifically, our theory is 
presented in three levels of complication; (i) the semiclas- 
sical theory of the Abelian Berry connection/curvature, 
which can be applied to 2D TIs in the adiabatic limit 
(Sec. |mAl ), (ii) the requantized effective theory of the 
general non-Abelian Berry connection/curvature, which 
can be applied to 2D as well as 3D TIs in the adiabatic 
limit (Sec. [17171] ), and (iii) the full quantum theory us¬ 
ing the Floquet Green’s function formalism, which can 
be applied to general situations (Sec. Ill C| ). We present 
computational results in Sec. |IV[ proving that the wind¬ 
ing number of the WSL provides a direct manifestation 
of topological invariants in both 2D and 3D TIs. We 
conclude in Sec. [V| where we discuss the experimental 
feasibility of observing the winding number of the WSL. 


II. BAND TOPOLOGY 
A. Chern numbers in two dimensions 

Let us begin with 2D TIs, where perpendicular spin 
components are conserved. The Ham iltonian for 2D 
TIs has the following generic structure!^ 33 -^ H = 
Ev=t,^L^(k)fe with H^k) = ff f *(-k) and 

H,(k)= eA+ d k ..q^ %k i- j, a, 

where V’L = ( c La> c Lv) with c La bein S the electron 
creation operator with momentum k, spin cr, and orbital 
a. The physical meaning of a depends on the specific 
model. Concretely, a denotes the sublattice indices, A 
or F>, in the Kane-Mele (KM) model, and the conduc¬ 
tion/valence band indices, E\ or H\, in the Bernevig- 
Hughes-Zhang (BHZ) model. I n is the n x n identity ma¬ 
trix. a = ( a x ,a y ,(J z ) consists of the Pauli matrices. The 
dk vector has three components dk = (dk,^, dk, y , dk, 2 ), 
where the first two components can be combined as 
dk,± = dk,^ ± id^y. The band dispersion is given by 
£±(k) = ek =b |dk| and ek =b |d_k| for spin up and down, 
respectively, indicating that the system becomes insulat¬ 
ing when |dk| ^ 0. 

The 2D band topology is characterized by the total 
flux of the Berry curvature piercing through the entire 
Brillouin zone (BZ) for each spin component. Called the 
(first) Chern number, the total flux of the Berry curva¬ 
ture is equivalent to the wrapping number of the normal¬ 
ized dk field around the unit sphere, 

= = ~a [ dk x dk y dk • (^/c^dk x <9/^dk), (2) 

47r J BZ 

where dk = dk/|dkp^. In fact, the 2D band topol¬ 
ogy can be fully determined by examining the low- 
energy behavior of dk=K+q around K = 0 or other 
low-energy momenta, where dje+q can be generally ex¬ 
panded as (Aq x , ±Aq y , M + F>(g 2 + g 2 )). When a sin¬ 
gle low-energy point exists at K = G 3 -^, the band topol¬ 
ogy becomes non-trivial if M/B < 0 and trivial other¬ 
wise. In the Kane-Mele mode? 2 34 defined on the hon¬ 
eycomb lattice, there are two low-energy Dirac points 
at K = K ± = (|j, both of which should satisfy 

the non-triviality condition in order for the whole valence 
band to become topologically non-trivial. 

In the presence of the time-reversal symmetry, the 
spin-dependent Chern numbers are always opposite be¬ 
tween different spin components, which means that the 
2D band topology is fully characterized by the Chern 
number difference. Motivated by the analogy between 
the charge and time-reversal polarization (TRP), the 
Chern number difference can be alternatively computed 
in a discrete form, which is formulated in terms of the 
parities of the time-reversal operator, ^(= 1 , 2 , 3 , 4 )? at four 
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time-reversal invariant momenta (TRIM^P^: 

4 

(-ir D «n*> ( 3 ) 

i= 1 

where the 2D Z 2 invariant, z/ 2 D, is identical to the half 
of the Chern number difference computed in an integral 
form in Eq. (§: 

v 2D = ~ 2 — (mod 2). (4) 

As shown in the following section, the fact that the 2D 
topological invariant can be computed in a discrete form 
has played an important role in defining the 3D topolog¬ 
ical invariants. 


B. Z 2 invariants in three dimensions 

In three dimensions, different spin components are in 
general mixed. The Hamiltonian for 3D TIs has the fol¬ 
lowing generic structurd^® H = -^"(k^k with 

H( k) = e k I 4 + d k • T 

/ ^k — dk,3 dk, 4 0 dk,- \ 

_ dk,4 + dk,3 dk,- 0 

0 dk,-i- ^k— dk,3 — dk,4 ’ 

V dk,+ 0 —dk,4 ek + dk,3/ 

(5) 

where V 4 = (4a 1 > < L 2 > < L 3 > c L 4 ) with C L bein S the 
electron creation operator with momentum k on gener¬ 
alized orbital <a, which includes both spin and orbital de¬ 
grees of freedom. As before, the physical meaning of a de¬ 
pends on the specific model. In the model for BiSe-family 
materials, (oq, <a 2 , CK 3 , oq) = (P1+P2~ P1+1, P2~ l 

). I n is the nxn identity matrix. T = (Ti,^,^,^,^) 
consists of the Gamma matrices satisfying the Clifford 
algebra {T^Tj} = 25 ^- 14 . Since there are five Gamma 
matrices, the corresponding dk vector has also five com¬ 
ponents; d k = (d k ,i,dk, 2 ,dk, 3 ,dk, 4 ,d k , 5 )- Specifically, 
for BiSe-family materials, dk can be expanded around 
the T point as dk,± = dk,i ±idk ,2 — A\(k x dk ,3 — 

M + B 1 (kl + ky) + B 2 kl , d k ,4 — Mk z , and d k ,5 — 0. Simi¬ 
larly, 6 k can be expanded as 6 k — C-\-Di(kl~\-k y )-\-D 2 k^. 

In the presence of mixing between different spin com¬ 
ponents, the spin-dependent Chern numbers cannot be 
defined properly in 3D TIs. In fact, mathematically, 
the Chern number cannot be defined at all in three di¬ 
mensions. Fortunately, if both inversion and the time- 
reversal symmetries are present, the 3D band topology 
can be characterized by four Z 2 invariants, (z/q; ^ 1 , ^ 2 , ^ 3 ), 
which depend on the parities o f the time-reversal oper¬ 
ator, ch(=i,... ,8), at eight TRIIVP^^. Playing the most 
important role by governing the robustness of a given 3D 
TI, the strong Z 2 invariant, Vo, is defined as 
8 

(-ir = rP = (- 1 ) i ' 2D (- 1 ) i ' 2D '’ ( fi ) 
2 = 1 


where z/ 2 d and z/ 2 d' are the 2D Z 2 invariants of the 
inversion-symmetric 2D subspaces containing one set of 
four TRIM and the other, respective!^ 29 . It is worth 
mentioning that is proportional to the integral of the 
Chern-Simons three form over the 3D BZ, which is only 
quantized in the presence of both inversion and the time- 
reversal symmetrical. 

Equation (|6| implies that the 3D TI becomes a strong 
TI if and only if the band topology of the inversion- 
symmetric 2D subspace containing one set of four TRIM 
is opposite to that of the 2D subspace containing the 
other. In the cubic lattice, this means that the band 
topology of the k z a = 0 subspace should be topologically 
non-trivial if that of the k z a = it subspace is topologically 
trivial, and vice versa. Of course, this statement should 
be true regardless of whether we choose the k x a = 0 and 
7 r subspaces (or the k y a = 0 and it subspaces) instead 
of the k z a = 0 and it counterparts. In this example, the 
other three Z 2 invariants simply correspond to the 2D 
Z 2 invariants of the k x = 0, k y = 0, and k z = 0 sub¬ 
spaces: that is to say, v x = ^ 2 D(/c x =o), ^2 = ^ 2 D(/c y =o), 
and v 3 = v 2D ( kz =o)- 


III. WINDING NUMBER OF THE 
WANNIER-STARK LADDER 

A. Semiclassical theory of the Abelian Berry 
connection / curvature 

To appreciate how the energy spectrum of electrons 
under an electric field can reveal the band topology, it is 
instructive to first consider the semiclassical dynamics of 
an electron wave packet moving in the lattice. Here, for 
simplicity, we assume that the electron wave packet is en¬ 
tirely composed of the plane waves consisting in a single 
nondegenerate band, which is well separated from other 
bands in the full energy spectrum. Also, for the time 
being, let us focus on the 2D TI, where perpendicular 
spin components are conserved so that the Hamiltonian 
for each spin component is decoupled. In this situation, 
the Berry connection/curvature becomes Abelian. An 
extension to the general case of the non-Abelian Berry 
connection/curvature is discussed in Sec. 

Under an electric field, the electron wave packet in the 
lattice performs the Bloch oscillation^, whose dynamics 
is described by the semiclassical Lagrangian 39 

C( r, r, k, k) = ftk • r + hA n (k) • k - H n (r, k), (7) 

where r is the center position, hk is the mean crystal 
momentum, and l~L n { r ,k) = £ n (k) + eE • r is the semi¬ 
classical Hamiltonian for the n-th energy band. A n ( k) = 
(0 n (k)|iVk|0n(k)) is the Berry connection with 0 n (k) 
being the periodic part of the Bloch wave function for the 
n-th energy band. In the case of the two-band model in 
Eq. 0 . there are two degenerate energy bands for each 
spin component, £±, a ( k), and the corresponding Berry 
connections, A±, a ( k). From this forward, we focus on the 
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FIG. 1. Schematic diagrams illustrating the fundamental connection between the 2 D topological invariant and the winding 
number of the WSL. Figure[l](a) depicts the semiclassical situation, where an electron wave packet performs the Bloch oscillation 
in the 2D momentum space under an electric field E. In the presence of the Abelian Berry curvature 23k, the Bloch oscillation 
is affected by the anomalous velocity eE x B^/h. The semiclassical path of the Bloch oscillation (encircling the torus along the 
k\\ direction) is quantized via the Bohr-Sommerfeld quantization rule with an additional geometrical factor e* 7Zak( ^^, where 
7zak(&_0 is the Zak phase. Figure [l] (b) shows that the center positions of the WSL eigenstates (i.e., the center positions 
of the envelope wave functions plotted via green and orange dashed lines) are shifted by 7z a k(&±)a||/27r, which generates a 
corresponding shift in the energy spectrum of the WSL eigenstates. For clarity, only the wave function for the central WSL 
branch is shown. Figure |T] (c) illustrates that the trivial/non-trivial Chern number is directly manifested in the trivial/non¬ 
trivial winding number of the WSL as a function of k±. Figure |T] (d) shows that, with the time-reversal symmetry dictating 
that the spin-dependent Chern numbers are opposite between different spin components, i.e., C^ = —C 4 ., two separate sets of 
the spin-dependent WSL branches (distinguished by red and blue lines) wind oppositely in the topologically non-trivial phase, 
which is characterized by the 2D Z 2 invariant = (Cf — C 4 ,)/2 = 1 (mod 2). In the topologically trivial phase characterized 
by z^ 2 D = 0, there is no winding of the WSL. Note that the guide lines for the WSL branches are obtained from the actual 
solution of the semiclassical theory for the Bernevig-Hughes-Zhang model via Eq. (11). 


lower occupied band with spin up so that we simplify the 
notation by setting H(r, k) = H_ 5 ^(r,k), £k = £_^(k), 
and Aik = Al- ? ^(k). Note that the same analyses pre¬ 
sented below can be repeated for spin down. 

By applying the variational method to the above La- 
grangian, one can derive two equations of motion. First, 
h k = — eE, which tells us that the crystal momentum 
parallel to the electric field, hk \\, changes linearly in time, 
while the perpendicular crystal momentum, hk ±, remains 
fixed. Here, the charge of electron is defined to be — e. 
Second, hr = Vk£k + eE x 23k, where the former term of 
the right-hand side is the usual group velocity and the lat¬ 
ter is the anomalous velocity due to the (Abelian) Berry 
curvature 3^ = Vk x Aik- Without the anomalous veloc¬ 
ity, these equations of motion describe the usual Bloch os¬ 
cillation, which can be understood in terms of the Bragg 
scattering of a wave packet at the BZ boundaries. 

The anomalous velocity generates not only a bending 
of the Bloch-oscillation orbit, but also a shift of its cen¬ 
ter. Considering that the quantized mode of the Bloch- 


oscillation orbit is nothing but the WSL eigenstatd^El, 
this means that the WSL eigenstate centers are shifted 
by the Berry curvature effect. To see how this is possi¬ 
ble, let us subtract a total derivative h-^(k • r) from the 
Lagrangian in Eq. 0 , which generates the following new 
Lagrangian 

£'(r, k, k) = —HTZ( r, k) • k - H( r, k), (8) 

where 7£(r, k) = r — Aik can be interpreted as the canon¬ 
ical position conjugate to k. In this interpretation, 
the Bloch-oscillation orbit is quantized according to the 
Bohr-Sommerfeld quantization rule: 

j) dky • 7£(r, k) = 2irn (n G Z), (9) 

where C denotes a closed orbit with the constant energy, 
where k\\a\\ sweeps through the entire BZ between — tt 
and 7 r with k±a± fixed. Here, a\\ is the lattice constant 
of the projected unit cell along the parallel direction to 
the electric field. Meanwhile, a± is defined so that a\\a±_ 
is the area of the unit cell. 
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As a consequence of the Bohr-Sommerfeld quantization 
rule, the center position, or the polarization of the WSL 
eigenstates is quantized according to § c dky • r = 2irn + 
7zak (fc_i_), where the non-integer shift is called the Zak 

phas^El. 


work is that the proper topological invariant character¬ 
izing the entire 2D band topology is not the Zak phase 
itself, but rather the change of the Zak phase across the 
one-dimensional BZ of fcj_, i.e., the winding number of 
the WSL. 


7Zak(&±) = j> 


dk\\ • A, 


( 10 ) 


which is generally a function of k±. See Appendix |A| for 
computational details of 7z a k(&_L)- Now, by equating the 
energy of the WSL eigenstates with the averaged semi- 
classical Hamiltonian over C, § c dk\\H(Y, k), one can 
show that the energy spectrum of the WSL eigenstates 
is given by 


f„ WSL (fc ± ) = £(k ± ) +(n + eEa\\ 


( 11 ) 


where £(k±) = ^ § c dk\\£^. 

The Zak phase is related with the Chern number, i.e., 
the total flux of the (Abelian) Berry curvature piercing 
through the entire BZ, 


c = T [ d 2 k • Bk, 

^ J BZ 


( 12 ) 


which, by using the Stoke’s theorem, can be rearranged 
as followPl 


C = 


L 


71 ^ a± dk _l Qyzak (k±) _ Ayzak 


-ir/a± 


2i r 


dk ± 


2i r 


(13) 


Equation (13) indicates that, if C = ±1, the WSL index 


n goes to n± 1 after k±a± sweeps through the entire BZ. 
This means that the Chern number is equivalent to the 
winding number of the WSL across the BZ as a func¬ 
tion of k _|_. Figure [l] (a)-(c) provide schematic diagrams 
summarizing the discussions so far. 

With both spin components taken into account, the 
time-reversal symmetry dictates that C ^ = — Cp This 
means that two separate sets of the WSL branches for 
spin up and down should wind oppositely and cross each 
other at k_\_a± =0 and ±7r, which are TRIM. Moreover, 
following the similar logic predicting that the Kramers 
doublets should exchange partners in the helical edge 
stated^, one can predict that the WSL branches should 
also exchange their Kramers-doublet partners (i.e., the 
doubly degenerate WSL eigenstates at TRIM) in the 
topologically non-trivial phase, while not in the trivial 
phase. This prediction is fully confirmed, as shown in 
Fig. [l] (d) and Sec. |IV| 

Finally, it is worthwhile to mention that the Zak phase 
was previously used as a topological invariant to predict 
the existence of edge states based on the bulk-edge corre¬ 
spondence between the quantized value of the Zak phase 
and the existence of a localized edge statrf 44 . In 2D TIs, 
the Zak phase is not generally quantized except at TRIM 
(or equivalently inversion-symmetric momenta), where it 
becomes either 0 or 7r. One of the main points in our 


B. Requantized effective theory of the general 
non-Abelian Berry connection/curvature 

We now consider the general case, where the Berry 
connection/curvature is non-Abelian 45 . To this end, it is 
convenient to concentrate on the generic four-band model 
discussed in Eq. © for 3D TIs, which has four energy 
bands in total with the lower two being degenerate and 
separated from the upper two (also degenerate) bands by 
an energy gap. As before, for simplicity, we focus on the 
lower two degenerate energy bands, both of which are 
fully occupied at half filling. Generally, in the absence 
of conserved (pseudo)spin components, the Berry con¬ 
nection/curvature has off-diagonal terms between differ¬ 
ent energy bands. Mathematically, the Berry connection 
in general has the SU( 2 ) non-Abelian gauge structure 
with Mk,a /3 = (0-,a(k)|iVk|^-,^(k)), where </>_ ?a (k) and 
k) are the periodic parts of the Bloch wave func¬ 
tion in the lower two degenerate energy bands with a 
and P denoting the pseudospin indices, say, u and d. In 
the non-Abelian case, the Berry curvature is defined as 
£>k = Vk x Mk — iAk x Mk, where the second term does 
not in general vanish due to the non-commutative re¬ 
lationship between the Berry connections with different 
spatial coordinates. 

The semiclassical Lagrangian can be extended in the 
non-Abelian case as followP^kM 

C( r, r, k, k, 77,77) =ihrfrj + hk • r + h^A^r]) • k 

— W(r, k), ( 14 ) 

where r, k, and H(r, k) = £(k) + eE • r are all defined 
the same as before in Eq. (pi. A new addition is the 
77 = {Vu,Vd) T variables, which denote the band decom¬ 
position of the wave packet satisfying the normalization 
condition, ^2 a=u d \da\ 2 = 1- The equation of motion for 
r is obtained similar to the Abelian case with a modi¬ 
fication that the Berry curvature £>k is now replaced by 
the 77 -averaged non-Abelian Berry curvature 77’*"#kT7. The 
dynamics of 77 is governed by ihr] = eE • A^r]. The equa¬ 
tion of motion for k is the same as before in the Abelian 
case. 

The Bohr-Sommerfeld quantization rule is no longer 
applicable in the non-Abelian case. Consequently, it is 
not obvious how to quantize the Bloch oscillation and 
obtain the energy spectrum of the WSL eigenstates. To 
overcome this ob stacl e, we requantize the semiclassical 
theory as follow^^. First, we note that the canon¬ 
ical position 7 £(r, k) in Eq. can be generalized to 
be 7 £(r, k, 77) = r — rj^A^rj. Then, we promote 7 Z to 
a quantum variable, 7^, which is conjugate to k. Next, 
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we formally eliminate the 77 variables by promoting the 
Hilbert space of the requantized effective Hamiltonian to 
be defined in the 77 pseudospinor space, I 77 ). Finally, by 
replacing r with 1 T 2 = 7ZI 2 + *4 k in H(r,k), we obtain 
the requantized effective Hamiltonian 

T~LeS (k) = S\J.2 + eE • (m 2 + M k ), (15) 


by imposing the periodic boundary condition: rj( k^, k\\ + 
27r/a||) = ?7(kj_, fey), which leads to the following quanti¬ 
zation condition: 

1 r 27T / a \\ 

— J dk\\ ^ WSL (kj_) - £ k , - eEAw,\i\ = 2n7r > 

( 21 ) 


where = iV k since it is conjugate to k. 

The energy spectrum of the WSL eigenstates can be 
obtained by solving the eigenvalue equation of / H e ff(k): 

^eff(k)|r 7 ) = £ WSL (kj_)|? 7 ), (16) 

where the WSL eigenenergy, £ WSL (k^), is shown to be a 
function of the perpendicular momentum to the electric 
field, k^. Note that the parallel momentum fcy is not a 
good quantum number. Now, it is convenient to utilize 
the periodicity in the momentum space and perform the 
Fourier transformation with respect to fcy, which leads to 
the Floquet-type representation of the eigenvalue equa¬ 
tion: 

E [M%?(M±)+neEa n 6 aP 6 nm ] V?m = f WSL (kj> a „, 

f3m 

(17) 

where 

2t r 

M^( k x ) = U j (£ k <5 a/3 + eE • .4 k>a/3 ) 

= f” m (k ± )^+eE-^(k ± ) (18) 

with a, /? E {rq d} and n, m E Z. In general, this Floquet- 
type eigenvalue equation is solved via numerical digonal- 
ization. 

Since actual solutions are obtained numerically, it 
is not easy to figure out the precise dependence of 
£ WSL (k ± ) on k^. It is, however, possible to infer 
a general structure of £ WSL (kj_) based on the follow¬ 
ing two observations. First, if £ WSL (k_L) is a solution, 
£ WSL (k±) + neEa\\ with n being any integer is also a 
solution. This means that £ WSL (k^) has a ladder-like 
structure even in the non-Abelian case. Second, in the 
Abelian limit, i.e., when A^ a p = M k (Wg, £ WSL (k±) re¬ 
covers the semiclassical solution in Eq. ( pTj ). To see this, 
let us rewrite the eigenvalue equation in Eq. ( fl 6 | ) in a 
differential equation form: 

£k+eE (^ +A ’") 

where M k) || is the parallel component of M k along the 
electric field. The above equation can be solved formally 
by 

r?(k) = e -A^" «*fc||[£ WSL ( k x)-£ k »-eB>t k ,, l ] } ( 20 ) 

where k = (kj_,fc||) and k' = (kx,fc||). Equation p0| 
is not yet a complete solution since the energy eigen¬ 
value £ WSL (k±) is unknown. £ WSL (k_\_) is determined 


V (k) = £ WSL (kx)7?(k), (19) 


which can be rearranged as the semiclassical solution in 
Eq. Based on these two observations, we deduce a 

general form of £ WSL (kx) as follows: 

C SL (kx) = £(kx) + (n + eEa h (22) 


where it is assumed that both £ and 7 do not depend on 
E at sufficiently weak E, i.e., in the adiabatic limit. This 
means actually that £ = £ § c dk\\£\t. Meanwhile, 

7 reduces to 7z a k only if the Berry connection/curvature 
becomes Abelian. It is confirmed in Sec. II V Bl that the ex¬ 
act eigenvalue solution of the requantized effective Hamil¬ 
tonian is indeed well described by Eq. (22) (and is also en¬ 
tirely consistent with the full quantum theory presented 
in the following section). 

So far, we have explained how to obtain the energy 
spectrum of the WSL eigenstates in the general case 
of the non-Abelian connection/curvature. With all the 
complications due to the non-Abelian structure, however, 
it is unclear at this stage how the energy spectrum of the 
WSL eigenstates can manifest the band topology in three 
dimensions. Fortunately, in the presence of both inver¬ 
sion and time-reversal symmetries, the 3D band topology 
is characterized by the band topologies of the special 2D 
subspaces within the 3D BZ, where the Berry connec¬ 
tions/curvatures become Abelian. These special 2 D sub¬ 
spaces are those with the inversion symmetry, where the 
inversion parity plays the role of a good quantum num¬ 
ber guaranteeing the existence of a diagonalized basis. In 
other words, within these special 2D subspaces, the Berry 
connection satisfies the following Abelian condition, 


K,4] = o, 


(23) 


with i and j denoting two orthogonal directions in the 
2D subspaces. In the cubic lattice, these 2D subspaces 
are the ki~kj planes with k^a = 0 and 7 r, where (i,j, k) = 
(x, 7 /, 2 ), ( 7 /, z, x), or (z,x,y). Under this condition, one 
can always find an appropriate set of bases, via which 
both M k and A J k are diagonalized simultaneously. Then, 
the WSL eigenenergy, f^ SL (k^), for each component of 
the diagonalized basis, say, pseudospin < 7 , simply reduces 
to the Abelian version with a well-defined Chern num¬ 
ber. Consequently, the WSL has a well-defined winding 
number for each diagonalized pseudospin component in 
these 2D subspaces. 

Now, with the time-reversal symmetry dictating that 
the total Chern number is zero, the three weak Z 2 in¬ 
variants are simply the 2D Z 2 invariants of three 2D sub¬ 
spaces at k x = 0 , k y = 0 , and k z = 0 , which are half 
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General 2D subspaces: 
Non-Abelian Berry 
connection/curvature 


Inversion-symmetric 
2D subspaces 
containing TRIM: 
Abelian Berry 
connection/ curvature 


FIG. 2. Schematic diagram illustrating how the 3D topological invariants are directly manifested in the winding numbers of 
the WSL within the inversion-symmetric 2D subspaces (purple-colored planes) containing TRIM (denoted by red circles), where 
the Berry connections/curvatures become Abelian. Here, the 2D subspaces are chosen to be parallel to the k x -k y plane within 
the 3D BZ without loss of generality. According to Eq. the strong Z 2 invariant is determined by the two Z 2 invariants of 
the 2D subspaces containing one set of four TRIM at k z a = 0 and the other at k z a = 7r. In general 2D subspaces (gold-colored 
planes), the Berry connections/curvatures are non-Abelian and thus the WSL does not have well-defined winding numbers. 


the Chern number differences between diagonalized pseu¬ 
dospin components within the corresponding 2D sub¬ 
spaces. Similarly, according to Eq. (J6|, the strong Z 2 
invariant is determined by the two Z 2 invariants of 2D 
subspaces, for example, at k z a = 0 and 7r. As mentioned 
previously, the same strong Z 2 invariant is obtained re¬ 
gardless of whether one chooses the k x a = 0 and 7 r, the 
k y a = 0 and 7r , or the k z a = 0 and 7 r planes. 

In addition to the general argument for the Abelian 
condition above, one can explicitly show that the Abelian 
condition is precisely satisfied in the generic four-band 
model in Eq. © for 3D TIs. Actual calculations are a 
little bit messy, but the gist of why the Abelian condition 
is satisfied in this model can be revealed by examining 
the long-wavelength expansion of A\^ around TRIM. For 
convenience, let us concentrate on one of the TRIM at 
k = 0, where A\t is expanded as 

Ak - £k o z z x k + x k, (24) 


where the concrete forms of the coefficients, and 
are not important for the current purpose. With help 
of Eq. (24), it is now straightforward to show that the 
Abelian condition is satisfied for all three 2D subspaces 
of the k x -k y , k y -k z , and k z -k x planes around k = 0. 
Similar arguments can be given for other TRIM. As men¬ 
tioned, it can be explicitly shown without any expansion 
that the Abelian condition is precisely satisfied for the 
entire subspaces of inversion-symmetric planes. See Ap¬ 
pendix [B] for the proof of the Abelian condition in the 
generic four-band model. 

To summarize, the 3D topological invariants are di¬ 
rectly manifested in the winding numbers of the WSL 
within the inversion-symmetric 2D subspaces containing 
TRIM, where the Berry connect ions/curvatures become 
Abelian. In general 2D subspaces, the WSL does not 


have well-defined winding numbers. Regardless of the 
existence of well-defined winding numbers, however, it 
is always possible to predict the energy spectrum of the 
WSL eigenstates precisely by solving the eigenvalue equa¬ 
tion of the requantized effective Hamiltonian in Eq. ( p~6| ) . 
This means that the requantized effective theory by it¬ 
self can be applied to any situations including 2D/3D 
topological insulators even with time-reversal symmetry 
breaking terms. See Fig. [2] for a schematic diagram illus¬ 
trating the discussions so far in this section. 

In our theory, the existence of well-defined 3D topo¬ 
logical invariants depends crucially on the fact that the 
the Berry connections/curvatures become Abelian and 
thus the winding numbers of the WSL are well defined 
within the inversion-symmetric 2D subspaces containing 
TRIM. It is important to note that, similar to our the¬ 
ory, the logic behind dimensional increase for the deriva¬ 
tion of 3D Z 2 invariants is also based on the fact that 
the 2D topological invariants are well defined within the 
inversion-symmetric 2D subspaces containing TRIM. In¬ 
terestingly, the discrete representation of the strong Z 2 
invariant in Eq. ^ becomes equivalent to an integral 
representation formulated in terms of the Chern-Simons 
three form if both inversion and time-reversal symme¬ 
tries are present 37 . Concretely, the integral of the Chern- 
Simons three form over the 3D BZ is given as follows 23 : 


= -! 
8tt J 


6= — d 3 k e ijk Tr 




, (25) 


where the non-Abelian Berry curvature is given by B £ = 
| e ljk B]i, where B^ = d^A^ - d^A^ - i[A^,A^\ with 
i,jA G {x,y,z}. In general, 9 in Eq. ( p5| is not quan- 
tized in contrast to the Chern number. Fortunately, 
however, in a close analogy with the fact that the one¬ 
dimensional Zak phase is quantized in the presence of 










the inversion symmetr}®!, it can be proved that 0 is 
quantized and related with the strong Z 2 invariant uq 
in the presence of both inversion and time-reversal sym¬ 
metries 37 : 


vq = — (mod 2). 

7T 


(26) 


In a sense, 6 can be regarded as a 3D analog of the one¬ 
dimensional Zak phase. It is important to note that the 
presence of both inversion and time-reversal symmetries 
is crucial to guarantee the existence of a global basis of 
wave functions, which is necessary to prove Eq. (26). 


Finally, it is worthwhile to mention that the Z 2 in¬ 
variants can be also obtained via the Wilson loop , which 
is the non-Abelian generalization of the Zak phase fac- 


W = V exp ( i j) dky • Ak 


(27) 


where V is the path ordering operator, A^ is the non- 
Abelian Berry connection, and C denotes a closed path 
traced by ky. (Note that a different sign convention is 
used for the phase in Eq. (27) compared to RefsP^O) 
Concretely, the Z 2 invariants are related with the winding 
numbers of the phase of the eigenvalues of the Wilson 
loop across the one-dimensional BZ of k^ 50 . Considering 
that the phase of the eigenvalues of the Wilson loop is 
in turn related with the center position of the Wannier 
function, this method is highly reminiscent of ours, where 
the Z 2 invariants are related with the winding numbers 
of the WSL. Below, we explain exactly how these two 
methods are connected. 

To this end, it is important to understand physically 
what the Wilson loop means. The physical meaning 
of the Wilson loop was elucidated by Grusdt et al“, 
who have shown that the Wilson loop is nothing but the 
propagator describing the Bloch oscillation in the limit 
of strong electric field. To appreciate further what this 
means precisely, let us examine the propagator at general 
strengths of electric field (for a closed path C ): 


IA — V exp 


1 <p dky 


• 4k + e. E Hk 


(28) 


where is the Hamiltonian in the absence of electric 
field. In the specific situation with two degenerate en¬ 
ergy bands described in this section, See Ap¬ 

pendix C in Reffor a detailed derivation of Eq. ([28). 


Now, it is important to note that, while expressed in 
a different gauge (namely, the time-dependent vector po¬ 
tential gauge via the Peierls substitution), U is actually 
equivalent to the propagator of our requantized effec¬ 
tive Hamiltonian H e ff in Eq. (15) (expressed in the time- 
independent scalar potential gauge). In this context, a 
rationale behind the usefulness of the Wilson loop is un¬ 
derstood as follow. Since the main information on band 
topology is embedded in Ak, one may ignore from 
the argument of the exponential in Eq. (28) if one is 


only interested in the characterization of the band topol¬ 
ogy, not the detailed coherent dynamics of the Bloch 
oscillation. Roughly, diagonalizing the Wilson loop is 
equivalent to diagonalizing H e ff while ignoring the band- 
dispersion part. Therefore, the winding numbers of the 
phase of the eigenvalues of the Wilson loop should carry 
essentially the same topological information as those of 
the WSL studied in this work. 

Rigorously, however, H \ c can be ignored only in the 
limit of strong electric field, i.e., Ue=oo = VV, where the 
adiabatic condition, which is necessary for the very va¬ 
lidity of the Berry connection/curvature, is completely 
violated. Therefore, the winding numbers of the phase 
of the eigenvalues of the Wilson loop may not be phys¬ 
ically observable due to the contradiction between the 
adiabatic condition and the strong electric-field limit. An 
improvement of the Wilson-loop scheme is to use directly 
the coherent dynamics of the Bloch oscillation, which is 
governed by Ue^oo- As mentioned in Sec. hi an interfero¬ 
metric method combining the coherent Bloch oscillation 
with the Ramsey interferometry was recently proposed to 
measure the Z 2 invariants in optical lattices 29 . It is em¬ 
phasized, however, that our method, where the winding 
numbers of the WSL are observed in spectroscopic mea¬ 
surements, can be applied to condensed matter systems, 
where the phase coherence is not guaranteed. 


C. Full quantum theory 

The semiclassical theory and the subsequent requan¬ 
tized effective theory are valid under the condition that 
the Bloch oscillation energy Ml = eEa\\ is sufficiently 
smaller than the interband energy difference. This con¬ 
dition is nothing but the adiabatic condition for the va¬ 
lidity of the Berry phase so that electrons can remain in a 
given band during the Bloch oscillation without making 
transitions to other bands. In other words, the electric 
field should be sufficiently weak so that the non-linear 
effects are negligible. On the other hand, the Bloch os¬ 
cillation should be sufficiently faster than the electron- 
impurity scattering rate so that electrons can complete 
a full cycle of the Bloch oscillation before scattered off. 
In what follows, we assess if these two conditions can be 
met simultaneously in the full quantum theory. 

The energy spectrum of the WSL eigenstates is re¬ 
vealed as a series of sharp peaks in the density of states 
(DOS), which can be obtained as the imaginary part 
of the retarded Green’s function, ImG£(o;). There are 
two gauge choices for the implementation of an elec¬ 
tric field; (i) the time-independent scalar potential gauge 
with 0 = — E x and (ii) the time-dependent vector poten¬ 
tial gauge with A = —cEt. In this work, we choose the 
time-dependent vector potential gauge, where the spatial 
translation symmetry is formally present. In the formal 
presence of the spatial translation symmetry, the mo¬ 
mentum parallel to the electric field, fey, is a conserved 
quantity. However, the DOS becomes gauge-invariant 
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and physically meaningful only if it is integrated over fey. 
The physically meaningful semi-local DOS is obtained by 
integrating out ImG^cj) with respect to fcy. The Green’s 
function in the frequency domain, G£(u;), is obtained 
from its counterpart in the time domain, G£(£, £'), via 
the Fourier transformation. In the time-dependent vec¬ 
tor potential gauge, G£(t, t') is computed via the Peierls 
shift, hk — >• ftk — eEt, of the Hamiltonian iif(k), i.e., 
H( k) —)> i^(k — eE t/h). It is interesting to note that the 
Peierls shift plays a role of incorporating one of the two 
semiclassical equations of motion mentioned previously 
in Sec. IIII A1 

Below, we explain how to compute the (retarded) 
Green’s function of the time-dependent Ham iltonian via 
the Floquet Green’s function formalism ! 5 ^ 52 i In the Flo- 
quet Green’s function formalism, the semi-local DOS is 
computed as follows: 

Pa (k ± ,u + nn) = -hmJ2(Gl) n M, (29) 

k\\ 

where the Floquet Green’s function is given by 

(Gk)S?M = J dt J dt , e i{uJ+nCl)t e- i{uJ+m ^ t ' (G k ) r a0 (t,t') 

(30) 


with —D /2 < uj < D/ 2 . Here, the Floquet frequency D is 
the natural frequency of the Peierls-shifted Hamiltonian, 
H(k — eE t/h), which depends on the specific structure 
of the Hamiltonian and is not necessarily the same as 
D. Usually, p a is summed over a since we are inter¬ 
ested in the semi-local DOS contributed by all general¬ 
ized orbitals including both spin and orbital degrees of 
freedom. However, in 2D TIs, where perpendicular spin 
components are conserved, p a for each spin component 
can be meaningful. Furthermore, in the special case of 
the Kane-Mele model defined on the honeycomb lattice, 
where the sublattice index is spatially distinguishable, p a 
for each spin and sublattice index can be independently 
meaningful. 

The retarded Green’s function in the time domain, 
(^k)a^(^,O, i s obtained by solving the following equa¬ 
tion: 

- l -Y,H a7 (k- em/h)(G r Mt,t'), 

7 

(31) 


where H a p( k — eE t/h) is the matrix element of the 
Peierls-shifted Hamiltonian between generalized orbital a 
and /3 at momentum k. Moving to the frequency domain 
via the Fourier transformation, Eq. (31) can be written 
in the Floquet matrix form: 


| + nD) + ir 7 

h 7 


Oa'y Onl 


Tjnl 


(k)}(G^M 




(32) 


where k) = ± / Q T dte^-^H^ik - eE t/h) with 

T = 27r/D. See Appendix |c| f or details on how to com¬ 
pute (Gk)™™ (<jo) from Eq. (32) for various TI models. 

The Floquet Green’s function formalism provides a 
natural platform to study the effects of electron-impurity 
scatterinjJ^. To this end, we take a simple model Hamil¬ 
tonian for the non-magnetic on-site electron-impurity in¬ 
teraction, 


H = V 




(33) 


where V is the electron-impurity interaction strength, 
and rii a and 7 ii mP; * are the electron and impurity number 
operators at the i-th lattice site, respectively. The full 
Green’s function, G£, is obtained by solving the Dyson 
equation 


( rnr—l\nm 

JaP 


r—l\nm 
a.p 


(^r 1 ) 


(^) n a n J a pS nm , 


(34) 


where (G£ )™™ is the inverse of the non-interacting 

Green’s function computed in Eq. (32) and E r is the self¬ 


energy due to the electron-impurity interaction. In this 
work, the impurity self-energy is computed via the self- 
consistent Born approximation (SCBA): 


(V)™M = ]T(gd™M, (35) 

k 

where Ui mp = y/n{ mp V with fi[ mp being the average im¬ 
purity number per site. It is important to note that, in 
principle, this formalism can be also applied to strongly 
correlated topological insulators, once the accurate self¬ 
energy is obtained for a strong electron-electron interac¬ 
tion. 

In what follows, it is shown that the results of the 
full quantum theory are entirely consistent with those 
of the semiclassical theory in the Abelian case and the 
requantized effective theory in the general non-Abelian 
case with all being robust against interband interference 
as well as non-magnetic impurity scattering. 


IV. RESULTS 
A. 2 D TI 

1. Bernevig-Hughes-Zhang (BHZ) model 

We study the BHZ modeP^as a first example of the 2D 
TI model with conserved perpendicular spin components. 
Imposed by the symmetry of the underlying microscopic 
structure around the V point, the Hamiltonian for four 
low-lying states (E\ t, Hi t? Hi |) can be written as 
the generic form in Eq. ([!]), which is repeated here for 

convenience: H = S k ,< 7 =t,iwith ffj,(k) = 
H*{- k) and 

H t (k) = + d k ■ = (7+7 et 7 t ,) • < 36 > 
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FIG. 3. Evolution of the semi-local DOS at (a) k±a = ±tt 
and (b) k±a = 0 as a function of electric field. Note that 
fan-shaped series of the WSL branches emerge from the cen¬ 
ter of the valence band, being entirely consistent with the 
semiclassical theory in Sec. |III A| Here, we consider a topo¬ 
logically non-trivial phase in the BHZ model with model pa¬ 
rameters such as A/AD — 0.6, B/AD = 0.6, C/AD — 0, and 
M/AD = —0.3. The red dashed lines are the semiclassical 
guide lines for the WSL branches obtained from Eq. HD- 

where VT = (4 b iCT with c L<t bein S the electron 

creation operator with momentum k and spin a (=t?i) 
on orbital a (= E\,H\). Around the T point, ek and 
dk can be expanded as ek — C + D{k 2 + k 2 ) and dk ~ 
(Ak x , — Ak y , M + B(k 2 + ky)), respectively. 

The low-energy Hamiltonian can be promoted to a 
tight-binding Hamiltonian via the minimal lattice reg¬ 
ularization, which replaces k by sin(A:a)/a and k 2 by 
2[1 — cos {ka)]/a 2 . Specifically, after the minimal lat¬ 
tice regularization, we set ek = C + 2D[2 — cos(k x a ) — 
cos (k y a)], dk,± = A[sin(A; x a) =F ism(k y a)], d^ z = M + 
2B[2 — cos (k x a) — cos (k y a)\, where we define A = A/a, 
B = B/a 2 , and D = D/a 2 , which all have the same phys¬ 
ical unit as M, i.e., energy. Here, we set the electric field 
to be aligned along the principal direction of the square 
lattice so that a\\ = a±_ = a. As mentioned previously, 
the band topology becomes non-trivial if M/B < 0 and 
trivial otherwise. 

Figure [3] shows the evolution of the semi-local DOS 
as a function of electric field, which exhibits fan-shaped 
series of the WSL branches emerging from the center 
of the valence band. As one can see, there is excellent 
agreement between the results of the semiclassical theory 
via Eq. for the Abelian Berry connection/curvature 
and those of the full quantum theory via Eq. ( [29] ) . It 
is important to note that, in addition to the main WSL 
branches emerging from the center, there are other WSL- 
like branches emerging from the band edges, which is 
known as the Franz-Keldysh effect!^! Similarly, the con¬ 
duction band (not shown in the figure) generates its own 
WSL eigenstate branches, which interfere with the WSL 


0 I ]8 



k±a± k±a± 

FIG. 4. Semi-local DOS as a function of k± in the BHZ 
model. Here, the semi-local DOS is summed over both spin 
and orbital degrees of freedom. The electric field is applied 
along the principal direction of the square lattice with mag¬ 
nitude eEa\\/AD — 0.16, which corresponds to the situations 
indicated by the red arrows in Fig. [3] Panels (a), (c), and 
(e) denote when the band topology is trivial with a choice of 
M/AD — 0.3, while panels (b), (d), and (f) denote when it 
is non-trivial with M/AD — —0.3. Here, A/AD, B/AD , and 
C/AD are chosen as the same as those in Fig. [ 3 ] The electron- 
impurity interaction strength is changed so that Vimp/4-D = 
0, 0.05, and 0.08 in the top [(a) and (b)], middle [(c) and 
(d)], and bottom [(e) and (f)] panels, respectively. The red 
and blue dashed lines denot e th e guide lines obtained in the 
semiclassical theory via Eq. HD for spin up and down, respec¬ 
tively. As predicted, the Kramers doublets exchange partners 
between k±a± = 0 and ±7r in the topologically non-trivial 
phase, while not in the trivial phase. 


eigenstate branches emerging from the valence band at 
sufficiently strong electric fields. Fortunately, despite all 
these complicated interferences, the main WSL eigen¬ 
state branches emerging from the center of the valence 
band can be clearly identified at an appropriate window 
of electric field, say, eEa\\/AD ~ 0.16, which is indicated 
by the red arrows in Fig. [3] 

Figure [ 4 ] shows the semi-local DOS at eEa\\/AD = 0.16 
as a function of k_\_, which confirms that the 2D topologi- 
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cal invariant is directly manifested in the winding number 
of the WSL, precisely as predicted by the semiclassical 
theory in Eq. (11). Specifically, two separate sets of the 
spin-dependent WSL branches wind non-trivially and op¬ 
positely in the topologically non-trivial phase (right pan¬ 
els) accompanied by an exchange of the Kramers-doublet 
partners between k±a± =0 and ±7r. Meanwhile, there is 
no winding of the WSL in the trivial phase (left panels). 

Now, to test the robustness of the band topology 
against non-magnetic impurity scattering, we investigate 
how the semi-local DOS changes as a function of electron- 
impurity interaction strength Ei mp . Note that, here, the 
effects of non-magnetic impurity scattering are taken into 
account within the SCBA via Eq. (35). As one can see 
from Fig. [4] (c)-(f), the winding number of the WSL can 
be clearly identified, unless the electron-impurity interac¬ 
tion strength becomes too strong to become comparable 
to the Bloch oscillation energy. 


2. Kane-Mele (KM) model 


0 I HI 8 



k±a± 


k±a± 


Next, we study the KM modelL^LLy, whose Hamiltonian 
is defined on the honeycomb lattice as follows: 

H = + i\ so E (37) 

(id) ((id)) 

where = (cj^, cj^), t is the hopping constant, and Vij = 
±1, depending on if the electron takes the left or right 
turn to get to the next-nearest neighbor. 

After the Fourier transformation, Eq. ( |37| ) reduces to 
the generic 2D TI Hamiltonian in Eq. (W) with e k = 
0, d k ,± = -t(e^ Cl + e^ ik - C2 + e^ ik - c % and d Kz = 
2 Aso[sin(k • ai) — sin(k • a 2 ) ~ sin(k • (ai — a 2 ))], where 
ci = a( 1/2, v%/2), c 2 = a(l/2,-v / 3/2), c 3 = a(-l,0), 
ai = a( 3/2, <\/3/2), and a 2 = a( 3/2,—\/3/2). As men¬ 
tioned previously, the band topology can be determined 
by examining the low-energy behaviors of d^ z around 
the Dirac points K ± = (|^, ± ^-); d kj * ~ t3\/3Aso ± 

fv^AsoK ^^) 2 + {q y a) 2 } at k = K ± + q. Low-energy 
behaviors of d^ z near both satisfy the non-triviality 
condition if Aso 7 ^ 0, which means that the whole va¬ 
lence band becomes topologically non-trivial under this 
condition. 

Figure [5] shows the comparison between the semi-local 
DOS of a topologically non-trivial phase (left panels) 
with Aso 7 ^ 0 and ordinary graphene (right panels) with 
Aso = 0, which again confirms that the 2D topological 
invariant is directly manifested in the winding number of 
the WSL. It is interesting to note that graphene does not 
have a well-defined value for the 2D topological invari¬ 
ant since the monopole singularities at K ± are located 
right within the 2D BZ. Figure [ 5 ] (b) and (d) show that 
the WSL energy spectrum behaves irregularly, when the 
semiclassical trajectory at a given k± passes through the 
monopole singularities at K^, making the Zak phase dis¬ 
continuous 44 . In some sense, graphene can be regarded 


FIG. 5. Semi-local DOS as a function of k± in the KM model. 
Here, the semi-local DOS is shown only for a given sublattice 
of the honeycomb lattice, while summed over the spin degree 
of freedom. The electric field with magnitude eEa\\/t = 0.16 
is applied along the armchair direction in panels (a) and (b), 
while along the zigzag direction in panels (c) and (d). See 
insets to see how the electric field is aligned in the BZ. Pan¬ 
els (a) and (c) correspond to a topologically non-trivial phase 
with Aso /t m 0.1, while panels (b) and (d) correspond to ordi¬ 
nary graphene with Aso = 0. The red and blue dashed lines 
denote the guide lines obtained in the semiclassical theory 
via Eq. (11) for spin up and down, respectively. The yellow 
dotted lines indicate when the semiclassical trajectory passes 
through the Dirac points, i.e., the monopole singularities. 


as being topologically critical, neither being topologically 
trivial nor non-trivial. 

Now, it is important to check that the winding number 
of the WSL does not depend on the electric-field direc¬ 
tion, while the detailed k ±-dependence of the WSL en¬ 
ergy spectrum may. To this end, in Fig. [ 5 ] (a) and (b), 
the electric field is applied along the armchair direction 
with a\\ = 3a/2 and a± = x/3a, while, in Fig. [H] (c) and 
(d), applied along the zigzag direction with ay = x/3a/2 
and a± = 3a. As one can see, the winding number of 
the WSL does not depend on the electric-field direction, 
being consistent with the fact that the winding number 
of the WSL is a topological quantity. 


B. 3D TI 

Finally, we study the 3D TI model describing strong 
3D TIs occurring in BiSe-family materials 3 4 . Around the 
T point, the Hamiltonian for four low-lying states (PI 4- 1 
, P 2“ t, PlJ P2“ 1) can be written as the generic form 
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in Eq. which is repeated here for convenience: H = 
^k^k^( k )^k with 

H( k) = e k I 4 + d k • T 

/ e k — d k? 3 d k , 4 0 
d k , 4 e k + d k , 3 d k5 _ 

0 d k , + e k - d k?3 

V d k>+ 0 —d k? 4 



(38) 


where Vl = ( 4 a 1 > 4 a 2 > 4 a 3 > c koJ with C L bein S the 
electron creation operator with momentum k on general¬ 
ized orbital a = (P1+ *f, P 2 “ *f, P 1 + P 2 “ |). As men¬ 

tioned before, d k can be expanded around the T point as 
^k,± = ^k,l A ^k,2 — i iky)) d k? 3 — M + L?1 (^T + 

ky) + P 2 ^, d k? 4 ~ A 2 & 2 , and d k ,5 ~ 0. Also, e k can be 
expanded similarly as e k cx C+D i(/^ + /^)+P 2 ^Simi¬ 
lar to before, the low-energy Hamiltonian in Eq. ( [38] ) can 
be promoted to a tight-binding Hamiltonian via the min¬ 
imal lattice regularization; e k = C + 2Ph[2 — cos (k x a) — 
cos (k y a)\ + 2 P 2[1 — cos(fc^a)], d k) ± = Ai[sin(fc x a) =b 
i sin (k y a)\ , d k ,3 = M + 2B\ [2 — cos (k x a) — cos (k y a)\ + 
2P 2 [1 — cos(A^a)], d k> 4 = A 2 sin (/^a), and d k? 5 = 0, 
where we define Ai = A^/a, Bi = Bi/a 2 , and Di = Di/a 2 
(z = 1,2). 

As mentioned previously, the 3D band topology is char¬ 
acterized by four different Z 2 invariants, (y 0 ; ^ 2 , ^ 3 )- 
Governing the robustness of a given 3D TI, the strong Z 2 
invariant vq becomes non-zero if the 2D topological in¬ 
variant of the inversion-symmetric 2D subspace contain¬ 
ing one set of four TRIM is different from that containing 
the other set. What this means in the cubic lattice is that 
the 2D topological invariant, or equivalently the winding 
number of the WSL, in the 2D subspace lying parallel to 
the k x -k y plane at k z a = 0 should be opposite to that 
at k z a = dz7r. Of course, this statement should be true 
regardless of whether we choose the k x a = 0 and it sub¬ 
spaces (or the k y a = 0 and n subspaces) instead of the 
k z a = 0 and it counterparts. 

The choice of the k z a = 0 and it subspaces is particu¬ 
larly convenient since, with = 0 at k z a = 0 and ± 7 r, 
the Hamiltonian in Eq. (38) becomes explicitly block- 
diagonalized with each 2x2 block precisely reducing to 


the 2D TI Hamiltonian of the BHZ model in Eq. (36). 
In this situation, exactly the same Zak phase analysis 
used for the 2D TI model can be applied to determine 
the winding number of the WSL in the 2 D subspaces at 
k z a = 0 and ±tt. Consequently, in the current model, 
the 3D band topology becomes non-trivial if the winding 
number of the WSL at k z a = 0 is non-zero while that 
at k z a = 7 T is zero, or vice versa. As one can see in 
Figs. [ 6 ] (a) and (f), this condition can be satisfied for an 
appropriate set of model parameters, generating a strong 
3D TI phase. 

Meanwhile, in general 2 D subspaces at k z a 7 ^ 0, ± 7 r, 
the Berry connections/curvatures become non-Abelian. 
In this situation, the Zak phase cannot be properly de¬ 
fined and thus the energy spectrum of the WSL eigen¬ 
states is no longer described by the simple semiclassi- 


0 I ]8 



-jt 0 jt -jt 0 jt 


k±a± k±a± 


FIG. 6. Semi-local DOS as a function of k± in the 
3D TI model within various 2D subspaces lying parallel to 
the k x -ky plane at different k z . Here, the semi-local DOS 
is summed over for all generalized orbitals including both 
spin and orbital degrees of freedom. The electric field with 
magnitude eEa\\/AD\ — 0.16 is applied along the princi¬ 
pal lattice direction within the x-y plane, in which situa¬ 
tion a\\ — a _l = a. Model parameters are chosen such that 
Ax/ADx = 0.6, A 2 /4Di = 0.5, Bi/ADx = 0.6, H 2 /TD 1 = 0.3, 
C/ADi = 0, D 2 /4Di = 0.2, and M/ADi — —0.3. Panels 
(a) and (f) describe the inversion-symmetric 2D subspaces, 
where the Berry connection/curvature becomes Abelian and 
thus the winding number of the WSL is well defined. The 
Kramers doublets exchange partners in panel (a) while not in 
panel (f), which means that the strong Z 2 invariant is non¬ 
trivial. In general 2D subspaces [panels (b)-(e)], the Berry 
connections/curvatures are non-Abelian. The energy spectra 
of the WSL eigenstates are accurately captured by the dashed 
guide lines obtained from the requantized effective theory via 
solving Eq. |l6| , which covers the Abelian situations in panels 
(a) and (f) as limiting cases. 


cal theory of the Abelian Berry connection/curvature 
in Sec. m Aj Instead, the energy spectrum of the 
WSL eigenstates is computed via the requantized ef¬ 
fective theory of the general non-Abelian Berry connec¬ 
tion/curvature in Sec. IIIB Figures [ 6 ] (b)-(e) show the 
comparison between the semi-local DOS obtained from 
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- Jt 0 JT -JT 

k±a± 


0 Jt 

k±a± 


FIG. 7. Counterpart of Fig. [6] within various 2D subspaces 
lying parallel to the k y -k z plane as a function of k x . Model 
parameters are the same as in Fig. [6] 


the full quantum theory and the guide lines of the WSL 
eigenstates (dashed lines) obtained from the requantized 
effective theory in Sec. Bj As one can see, the agree¬ 
ment between the two theories is excellent. It is im¬ 
portant to note that the requantized effective theory re¬ 
produces the results of the semiclassical theory of the 
Abelian Barry connection/curvature in Fig. [6] (a) and (f) 
in the Abelian limit. 

Finally, it is important to check if the strong Z 2 in¬ 
variant is uniquely determined, being independent of the 
choice of 2D subspaces. To this end, we examine the 
evolution of the semi-local DOS in various 2D subspaces 
lying parallel to the k y -k z plane as a function of fc x , which 
is shown in Fig. [7| Unlike in Figs. [6] (a) and (f), here, 
the Hamiltonian is not explicitly block-diagonalized even 
within the inversion-symmetric 2D subspaces at k x a = 0 
and 7 r. Fortunately, as explained in Sec. EH the Berry 
connections satisfy the Abelian condition in Eq. (23) 
within these 2D subspaces. This means that, within 
these 2D subspaces, the Hamiltonian is decomposed into 
two independent parts with each having the well-defined 
winding number of the WSL. Moreover, the time-reversal 
symmetry dictates that the two winding numbers should 


be opposite. As one can see, this is exactly confirmed in 
Figs. [7] (a) and (f). As before, in general 2D subspaces 
where the Berry connections/curvatures are non-Abelian, 
the energy spectrum of the WSL eigenstates is accurately 
captured by the requantized effective theory in Sec. EH 
A similar analysis for the 2D subspaces lying parallel to 
the k z -k x plane is guaranteed to generate exactly the 
same conclusion as the above due to the reflection sym¬ 
metry between the x and y directions. 


V. DISCUSSION 

In this work, we show that the non-trivial band topolo¬ 
gies of both 2D and 3D TIs, characterized by the Chern 
numbers and the Z 2 invariants, respectively, are directly 
manifested in the winding numbers of the WSL emerg¬ 
ing under an electric field. Being alternative to the topo¬ 
logical magneto-electric effecl 2 ^, this provides a spectro¬ 
scopic method to measure the topological invariants di¬ 
rectly in the bulk of both 2D and 3D TIs. Below, we 
discuss briefly how this method can be realized in actual 
experiments. 

The main physical observable to be measured is the 
semi-local DOS. Considering that the modern STM tech¬ 
nique has the sufficient spatial resolution to distinguish 
individual atoms in a crystal, the semi-local DOS can 
be in principle obtained by scanning the surface of a 2D 
TI (possibly, obtained in the thin-film limit of 3D Tld^J) 
or the cleaved surfaces of a 3D TI and then partially 
Fourier-transforming the STM data along the perpendic¬ 
ular direction to the electric field. Another method to 
measure the semi-local DOS is the ARPES, which can 
give rise to the momentum-resolved information directly 
near the cleaved surfaces of a 3D TI. 

The WSL has been so far observed only in man-made 
structures such as optical and semiconductor superlat¬ 
tices since the typical lattice constant in a natural crystal 
is usually too small (^ a few A) that the energy spectrum 
of the WSL eigenstates (broadened by impurity scatter¬ 
ing) is not well resolved for a typical strength of electric 
fiel d 40 * 41 1 . A key task is to apply a sufficiently strong 
electric field to overcome the broadening effect, while 
suppressing the Joule heating. This has been achieved 
in semiconductor superlattices with the superperiod of 
50-100 A at an electric field of 10-20 kV/cm, in which sit¬ 
uation the energy spacing between the WSL eigenstates 
is roughly 10-40 meV 40 41 4 If experiments can be per¬ 
formed in a natural crystal with the same energy resolu¬ 
tion, the required electric field is estimated to be roughly 
in the order of 100-200 kV/cm. With stronger electric 
fields, it would be important to make TIs truly bulk- 
insulating 55 57 in order to suppress the Joule heating. 

Actually, it has been proposed that an artificial TI (as 
well as an artificial Weyl semimetal) can be constructed 
in a superlattice structure, which is composed of alter¬ 
nating layers of topological and ordinary insulators with 
the layer thickness spanning many ( rsj 20-30) unit celld^. 
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In this situation, the energy spacing between the WSL 
eigenstates can be dramatically enlarged, opening up the 
possibility of observing the winding number of the WSL 
even in a weak electric field. For future work, we would 
like to investigate if this possibility can be realized in 
realistic material condition! 59 !. 

Furthermore, it would be also interesting to observe 
the energy spectrum of the WSL eigenstates in graphene, 
which is topologically critical as explained in Fig. [5] In 
graphene, the limitation of a small unit cell can be ef¬ 
fectively overcome by forming the moire structure®®^. 
Moreover, extensive efforts are being devoted to enhance 
the spin-orbit coupling strength ^ 3 to realize the KM 
model in actual graphene. 
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Appendix A: Zak phase 


In this section, we provide computational details of 
the Zak phase in 2D TIs, where the Berry connec¬ 
tion/curvature is Abelian. The 2D TI Hamiltonian 
in Eq. 0 is repeated here for convenience: H = 
with = H t *(-k) and 


fft(k) = 


_ ( dk,- 


d] 


k,+ 


e k dk,2 


(Al) 


where the detailed form of dk is not important for the 
current purpose except that |dk| = |d_k| due to the time- 
reversal symmetry. 

After diagonalizing Eq. (Al), the energy eigenvalue is 
obtained as follows: 


where Ck,± = (d k> * ± |d k |)/ \J (d k ,*) 2 + (d k , y ) 2 and <p k = 
tan (dk,y/dk,x) • 

Now, with help of Eqs. ( A3| ) and ( |A4| ), the Berry con¬ 
nections for spin up and down can be computed as fol¬ 
lows: 


At >t (k) = (0±,t(k)|iV k |0±,t(k)> 

= C^k,± (dk, 2 / Vkdk,cc d\^ x Vkdk,y) •, (A5) 
= (t±M |iVk|^±4(k)) = -At lt (k) (A 6 ) 


where 


k,± 


_I dk | A dk,z_ 

2|d k | [(dk,cc) 2 + (dk, y ) 2 ] 


(A7) 


As discussed in the main text, the energy spectrum of 
the WSL eigenstates for each spin component a is given 
by 


au k ±) =^±(^)+(«+ 7Zsk ’^ (fcx) ) eEa b 

J (A 8 ) 

where £±(k±) = § c dk\\£±(k). The spin-dependent 
Zak phase, 7 zak,±,cr (&±) 5 is evaluated via 


7Zak,±,<r(fc_L) = dk|| ' A± J<T (k), 


(A9) 


whic h in dicates that 7zak,±,t(fc-i_) = — 7zak,±,t(^_L) due to 
Eq. (A 6 ), which in turn means that the winding numbers 
are opposite for different spin components. 

It is important to note that the same formalism can 
be applied to the inversion-symmetric 2D subspaces 
within the 3D BZ of 3D TIs, where the Berry connec¬ 
tions/curvatures are Abelian. By choosing the right basis 
of wave fun ctions, the 3D TI Hamiltonian can be written 
as Eq. with conserved pseudospin components. 


Appendix B: Proof of the Abelian condition 

In this appendix, we provide the proof of the Abelian 
condition in Eq. ( [23] ) for the generic 3D TI model. To 
this end, let us rewrite the generic 3D TI Hamiltonian in 
Eq. (§: H = £ k <4#( k )</>k with 


£±(k) = £k ± |dk|, (A2) 


for both spins, indicating that the system becomes insu¬ 
lating when |dk| 7 ^ 0 in the entire BZ. The corresponding 
eigenstates for the lower and upper bands, </>_ )Cr (k) and 
0 + 5 Cr (k), respectively, are given as follows: 


k±,t( k )> 

I0±4(k)> 


1 ( Ck,±e _ * Vk 

yi + (Ck ,±) 2 v 1 

1 / —Ck,±e* Vk 

71 + (Ck .±) 2 V 1 


(A3) 

(A4) 


(£k — dk,3 dk,4 0 dk.- \ 

£k + dk,3 dk.— 0 

dk,+ £k - dk,3 —dk,4 
0 — dk,4 ek + dk,3/ 

(Bl) 

where, via the minimal lattice regularization of the low- 
energy effective model, the system parameters can be 
written as ek = C+2Di [2—cos (k x a)— cos (k y a)]-\-2D2 [1 — 
cos(fc^a)], dk,± = Ai[sin (k x a) ± i sin (k y a)], dk ,3 = 
M + 2Bi [2 — cos (k x a) — cos (k y a)\ + 2 B 2 [1 — cos (k z a )\, 
dk ,4 = A 2 sin (. k z a ), and dk,s = 0. 


H(k) = 


«k,4 

0 

^k,+ 
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Similar to the 2D TI case, the energy eigenvalue is 
given by 


£±(k) = e k ± |d k |, (B2) 


which indicates that there is a double degeneracy for both 
upper and lower bands. Let us distinguish the two de¬ 
generate energy eigenstates within the upper (subscript 
+) and lower (subscript —) bands by introducing a pseu¬ 
dospin index, say, u and d. In other words, the energy 
eigenstates are distinguished by two indices with one be¬ 
ing dz and the other being u/d. Specifically, the energy 
eigenstates, |</>±,u(k)) and |0± jd (k)), are given as follows: 


|0±,«(k)> 


yi + ( Xk ) 2 + (Ck,±) 2 


|0±,d(k)) = 


yi + (xk) 2 + (Ck,±) 2 


Xk e ~ tVk 

(k,± e ~ itPk 

1 
0 

1 
0 

Xk e^ k 

Ck,±e* Vk 


(B3) 

(B4) 


where Xk = <4,4/V(<4 ,i) 2 + (<4,2) 2 , Ck,± = (<4,3 ± 
|d k |)/V(dk,i) 2 + ( rf k, 2 ) 2 , and <p k = tan _1 (<4 >2 /<i k ,i). 

In general, the Berry connection has off-diagonal ma¬ 
trix elements mixing between |0±,u(k)) and |0± jd (k)), 
which generates the SU(2) non-Abelian gauge struc¬ 
ture^^. Specifically, the Berry connection for the lower 
band is explicitly written as follows: 


*^k,uu = (0_, w (k)|iV k |0_, u (k)) 

= - <a k (d k? 2 v k d k? 1 - d k; iV k d k?2 ), (B5) 

-4k, dd =(0-,d( k)|iV k |0_ 5 d(k)) 

=a k (d k ,2V k d k) i — d k? iV k d k)2 ), (B6) 

^k,ud =At 4 u = (0-, n (k)|iV k |0_,,(k)) 

=/3k [dk,4(~^V k d k? i — V k d k52 ) 

+ (idk,i + dk,2)V k d k? 4], (B7) 


where 

«k = 

A< = 


(|d k | — d k , 3 ) 2 + (d k , 4 ) 2 

2|d k |(|d k | - d k , 3 )[(<4,i) 2 + (d k)2 ) 2 ] ’ 1 ’ 

_ ~(<4,i ±idkrf _ mq x 

2|d k |(|d k |-d k , 3 )[(d k ,i) 2 + (d k ,2) 2 ]' 1 ’ 

6l), and (|B7|) as 


Now, one can rearrange Eqs. 
follows: 


-4 k =/3 k AiA 2 asin(£^a) cos (k x a)cr y 


— a^A 2 a sm{k y a) cos {k x a)cr z . 

(BIO) 

A^. = — /3 k AiA 2 asin(^a) cos (k y a)a x 


+ a^A 2 a sin(k x a) cos (k y a)cr z . 

(BH) 

-4 k =/3 k AiA 2 asin(/c ?/ a) cos (k z a)cr x 


— /3 k AiA 2 asin(/c x a) cos (k z a)cr y , 

(B12) 


where the Pauli matrices, (a x , a y: a z ) : are represented in 
the basis of |0_ jU (k)) and |0_ jd (k)) as follows: 


^ = |0_ >t4 (k))(0_, d (k)| + |0_, d (k))(0_, tt (k)|, (B13) 

a y = —i|</>-, u (k)) (</>-,d(k) | + i|<£-d(k))(<£_ u (k)|, 

(B14) 

Oz = l0-,ti(k))(0-,u(k)| - |^-,d(k))(^-,d(k)|. (B15) 


Then, after some algebra, one can show that the com¬ 
mutator between various components of the Berry con¬ 
nection is summarized compactly as follows: 


[*4 k , A^\ = Zitijk sin {k^a) cos {kid) cos (fc^-a) (B16) 


where 


Tk =a k/ d k A^A 2 a 2 [sin {k x a)cr x + sin (. k y a)a y \ 

+ /d^A^A^a 2 sin {k z a)cr z . (B17) 


Equation (B16) implies that the commutator vanishes in 
the entire ki~kj plane if = 0, ±7r/a, which is nothing 
but the Abelian condition for the inversion-symmetric 
2D subspaces. This completes the proof of the Abelian 
condition in Eq. (23). 


Appendix C: Floquet Green’s function 


The goal of this appendix is to explain how to compute 
the Floquet Green’s function, which satisfies the follow¬ 
ing equation: 


| ^h{u + nCl) + ir] 

1,1 


S ai 6 nl - H%( k)} (G r k ) l ^(u) 


= <W 5 


nmi 


(Cl) 


where H™™( k) = ^ J Q T dte l ( n 3 (k — eEt/H) with 

T = 27r/Cl. Since the concrete form of the Hamiltonian 
matrix element H a p( k) depends on the specific model, it 
is not possible to obtain the general solution for the Flo¬ 
quet Green’s function in a closed analytic form. Fortu¬ 
nately, considering the structure of the generic TI Hamil¬ 
tonian, it is possible to derive a formal solution for the 
Floquet Green’s function with the same orbital indices, 
(^k)So', by summing away all other contributions from 
those with different orbital indices. This formal solution 
is convenient since the semi-local DOS is solely dependent 
on the orbital-diagonal components (G k )^2- Below, we 
present such a formal solution first in the case of the 2D 
TI with conserved spin components and then in the case 
of 3D TI with mixed spin components. 
















16 


1. 2D TI with conserved spin components 


We begin by rewriting the 2D TI Hamiltonian in a 
matrix form including both spin components: 


H{ k) 


( dk,- 0 0 ^ 

dk,+ He — dk,£ 0 0 

0 0 6 —k T d— k ? £ d— k,+ 

V 0 0 d— k,— C—k d— k,^ / 

(C2) 


where the concrete forms of ek and dk = (d^ x , d^ y , d^ z ) 
depend on the specific model and are shown in Sec. IV A| 
As before, dk,± = d^ x ± id^ y . Here, the basis is chosen 
such that (1, 2, 3,4) = (E\ t, Hi t, E\ Hi |) for the BHZ 
model and (At, B\, A\,, B\) for the KM model. 

By summing away all contributions from the orbital- 
off-diagonal components in Eq. (Cl), the inverse of 
the orbital-diagonal components of the Floquet Green’s 
function can be written in a compact notation with 
[G Kaa ] nm = (Gk)2a as follows: 


Gk7n(w) = [P^M]" 1 - D- • Pj» • D+ (C3) 
GJ^M = [Pk H]" 1 - D+ • P+M • D k , (C4) 

GL; 3 3^) = [Pt k M] _1 - Dl k • Pl k (w) • Dl k , (C5) 
G£7 4 » = [P^M]" 1 - Dl k • P7 k H • D7 k , (C6) 

where the Floquet matrices P k (w) and D k are given as 


[p£M] 


— 1 ,nm 


[h(u + nfi) + irj\6 nm - (e£ m ± dl m z ), 

(C7) 


and 


( D ±)«m =d nm , ( C8 ) 

where e£ m = ± / Q T dte i{m - n ) ht e {k - eEt/h) and d£ m = 
d / 0 T dte*( m_n )^ t d(k - eEt/h) with T = 2 tt/0. 

The specific forms of e^ n and d™ depend on not only 
the detailed k dependences of ek and dk, but also the 
electric-field direction. In what follows, we present the 
specific forms of e^ n and d^ n for the BHZ model along 
the principal direction of the square lattice and the KM 
model along the armchair as well as the zigzag directions. 


a. BHZ model 

The BHZ model is defined on the square lattice. To 
appreciate that the concrete form of the Floquet matri¬ 
ces depends on the electric-field direction, let us imagine 
that the electric field is applied along the direction with 
angle 6 measured from the principal, say, x axis of the 
square lattice. In this situation, the Peierls-shifted crys¬ 


tal momentum becomes 


f hk x + ^A x {t) \ _ ( cos# — sin#\ f Kk\\ — eEt\ 
yhky + ^A y (t) J ~ l^sin# cos 0 J y hk± J 

_ 1 / (qhk\\ — phk±) — qeEt \ 

~ ^y p 2 + q 2 y (pfifcy + qhk±) - peEt J ’ 


where we set 0 = tan -1 (p/g) with p,q E Z so that 
all different time-dependent terms in e(k — eEt/h) and 
d(k — eEt/h) become commensurate with each other. 
In other words, the entire time dependence occurs 
through cos {k x a — eaA x (t)/hc ), cos (k y a — eaA y (t)/he ), 
and their sine counterparts, which means that there ex¬ 
ist two different oscillation frequencies: qeEa/y/q 2 + p 2 
and peEa /y/ q 2 + p 2 . In this situation, the Floquet fre¬ 
quency, which is the natural frequency of the Peierls- 
shifted Hamiltonian, is given as H = eEa\\/h with a\\/a = 

gcd(g, p)/\Jq 2 +p 2 , where gcd(g,p) denotes the greatest 
common divisor of q and p. Note that, for general angle 
6,0. is not necessarily the same as the Bloch oscillation 
frequency O = eEa\\/h. 

In the case of (q,p) = (1,0), where a\\ = ay = a and 
thus 0 = 0, the Floquet matrices are given as follows: 


£ nm = e i(n-m)k„a[ c + 2 £ (2 _ cos ( k±a ^ Snm 

~ -D(<5n,m+1 + <Vm-l)] , (CIO) 

d ran = ±_ e i(n-m)k n [ ± 2 S in(fc ± a)<W 

T ^n,m+l ,m—1 ], (cii) 

d mn = e i(n-m)k^a + 2B (2 - COs(k ± a))]5 nm 

~ + <5n,m-l)}, (C12) 


which can be plugged into Eqs. ( |C7| ) and (C8) to com¬ 
pute the i nver ses of the Floquet Green’s functions in 
Eqs (C3)-(C6), which are then inverted to generate the 
Floquet Green’s functions themselves. 


b. KM model 

The KM model is defined on the honeycomb lattice, 
where the Peierls-shifted crystal momentum is given by 

( hk x + ~A x (t) \ _ f cos 6 — sin#\ / hk\\ — eEt\ 
yhky + ^A y (t) J l^sin# cos# J y Hk± J 

= 1 f q{fdc\\ ~ eEt ) rri1 

a/(/ 2 + p 2 /3 \p{hk || - eEt)/V 3 + qk± J ’ 

where we set 9 = tan -1 (~^) p, q € Z so that 

(q,p) = (1,0) and (1,1) correspond to when the electric 
field is applied along the armchair and the zigzag direc¬ 
tions, respectively. Note that (q,p) = (0,1) also corre¬ 
sponds to the zigzag direction. The entire time depen¬ 
dence of the Peierls-shifted Hamiltonian occurs through 
six terms; exp (tk(t) • ci), exp (ik(t) • C 2 ), exp (ik(t) • C 3 ), 
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FIG. 8. Unit cells and the corresponding Brillouin zones 
in the electric-field-applied KM model, (a) The red and blue 
boxes denote the unit cells in the lattice when the electric 
held is applied along the armchair (red arrow) and the zigzag 
(blue arrow) direction, respectively. Here, A and B denote 
different sublattices, (b) The corresponding Brillouin zones 
are shown as the red and blue boxes for the armchair and 
the zigzag directions, respectively. Here, = (|^ , ± 3 ^ a ) 
denote the Dirac points. 


Second, the zigzag direction is obtained by choosing 
(g,p) = (0,1), where dy = ay = \/Zaj 2 and a± = 3 a. 
In this situation, the Floquet matrices are given by 


d nm e i(n-m)k\\“H[- exp (Tik±a ± /3)5 nm 

- exp (±ikj_aj_/ 6 )(S n>m+1 + i)], (C16) 

d nrn = S n - m ^ & nso[-2ic08{k ± aj2){S n , m+ l ~ 1 ) 

+ i{Sn,m +2 — 3n,m- 2 )]- (C17) 


Note that the unit cells for the armchair and the zigzag 
directions have different shapes, but the same area. See 
Fig.0 for illustration. 

2. 3D TI with mixed spin components 


sin (k(£) • ai), sin (k(£) • a 2 ), and sin (k(£) • (ai — a 2 )), 
where Ci = a( 1/2, \/3/2), c 2 = a( 1/2, — \/3/2), C 3 = 
a(— 1,0), 3.1 = a( 3/2, V3/2), and a 2 = a(3/2, — a/3/2). 
Note that k(£) = k — eEt/h. This means that there are 
six different oscillation frequencies. As a consequence, 
the Floquet frequency is given by Cl = eEa\\/h with 

a\\/a = gcd(3 q+p, 3 q-p,q+p,q-p, 2q, 2p)/2yY +p 2 /3, 
where gcd(ai, • • • ,a n ) denotes the greatest common di¬ 
visor among (ai,••• , a n ). 

In the KM model, is always zero. In this work, 
d^_ and d are computed in two different situations 
with the electric field applied along the armchair and 
the zigzag directions. First, the armchair direction is 
obtained by choosing (q,p) = (1,0), where dy = a/2, 
while ay = 3a/2 and a± = \/3 a. In this situation, the 
Floquet matrices d and d are given by 

d nm ± = [-2 cos (k ± a ± /2)5 n , mTl - S n , m±2 ] , 

(C14) 

d nrn = j{n-m )kf H Ago[ _ 2 sin ( k±a± )5 nm 

-(- 2 sin (k±a±/ 2 )( 6 ntm +3 -t- 5 n , m — 3 )]. (C15) 


Let us begin by rewriting the 3D TI Hamiltonian in a 
matrix form: 


H( k) 


/ e k ^k,3 ^k,4 0 ^k,— \ 

dk,4 Gc + dk,3 dk,- 0 I 

0 dk,+ ek — dk,3 —dk,4 I 

V dk,+ 0 -d k , 4 ek + d k , 3 / 

(CIS) 


where the basis is chosen such that (1,2,3,4) = (P1+ t 
, P 2“ t, Pl/~ Jr, P2“ i). The concrete forms of ek and dk 
are shown in Sec. hvh 

Similar to the 2D TI case, by summing away all con¬ 
tributions from the orbital-off-diagonal components in 
Eq. ( |C1| ), the inverse of the orbital-diagonal components 
of the Floquet Green’s function can be written as follows: 
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GL7i» = [p^MP* - D7 • QkH • D+ 

- (pi - D7 • Q£M • D 4 • PJ» • D+) • R£+M • (d 4 - ■ P~(co) • D k • Q£M • D+), (C19) 

GJ^M = [P^MD - D 4 • Q*+M • D 4 

- (d^ - D 4 • Q*+M • B~ • P+M • D 4 ) • r£» • (d+ - D 4 • P+M • D+ • Q»+M • D 4 k ), (C20) 

GJ^M = [P7(w)] _1 - D+ • QJ(w) • B~ 

- (d 4 - D+ • Q£(w) • D 4 • Pj» • D k ) • R£" M • (d 4 - D+ • Pj» • D 4 • Q£(w) • D k ), (C21) 
GJ^M = [p+MP 1 - D 4 • QpM • D 4 

- (d+ - D 4 • Qp (w) D+ P+ (w) • D 4 ) • Rp (w) • (d* - D 4 P+ (w) • D k • Qp (w) • D 4 ) , (C22) 


where 


[RpMP 1 

[RpMP 1 


[P+H]- 1 - DJ • P k (w) D+ - Dj Pj» • D 4 • QJ(w) • D 4 • P~(co) • D+ 
PPM” 1 - D 4 • P+M D 4 - D 4 P+M • D+ • QfH • D7 • P+M • Dk, 


(C23) 

(C24) 


[QkH ]" 1 = [p£MP* - D£ • p k M • Dk, (C25) 

[Qf H ]" 1 = [ P kM ] _1 - DJ • P+M • D+ (C26) 


where the Floquet matrices P k M and Dp 4 are given 
as 

[P+M]- 1 -"" 1 = [ft(w + nQ) + irf\5 nm - (e£ ro ± d£^), 

(C27) 

and 

(D+) nm =dl m ± , (C28) 

( D 4)nm =rf nm (C29) 


with (ek)mn and (dk) mn defined the same as before. 

As mentioned in the main text, we are interested in the 
winding number of the WSL within various 2D subspaces 
to determine the strong Z 2 invariant. As an example of 
such 2D subspaces, let us first consider the 2D subspaces 
lying parallel to the k x -k y plane with the electric field 
applied along the principal, say, x direction, in which 
case ay = ay- = a± = a. Then, the Floquet matrices are 
given as a function of two conserved momenta k±(= k y ) 
and k z as follows: 


e nm = e i{n-m)k,a^ C + 2 £> 2 ( 1 - C0S (k z d)) + 2^(2 - COS (k ± a))]S nm - + S n , m - 1 )}, (C30) 

dl m ± = [2 sin ( k ± a)6 nm T S n , m+ 1 ± S n , m - 1 ] , (C31) 

d nrn = jfr-m )*„«{ [ M + 2 B 2 (1 - COS (fc,a)) + 25,(2 - COS (k±a))]S nm ~ B 1 (S n ,m+ 1 + (C32) 

M Sin ( k z a)5 nm , (C33) 


which can b e plug ged into Eqs. (C27)—(C29) to compute 
Eqs. (C23)—(C26) and subsequently the inverses of the 
Floquet Green’s functions in Eqs. (C19)—(C22), which are 
then inverted to generate the Floquet Green’s functions. 
Similarly, we also consider the 2D subspaces lying par¬ 


allel to the k y -k z plane. Now, the electric field is applied 
along the y direction in these 2D subspaces. Then, the 
Floquet matrices are given as a function of two conserved 
momenta k_ l(= k z ) and k x as follows: 
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e nm = e i(n-m) k]] a { ^ + 2 £ 2 (1 _ cog (A;±a )) + 2^(2 - COS (fc*a))]<W - (£„, m+ l + (C34) 

d nm = l _ e i(n- m) k,a Ai [3 gin ± <y m+1 q= 6^^] , (C 35 ) 

d nrn = [ M + 2 S 2 (1 - C OS (k ± a)) + 2^(2 - COS (k x a))]S nm - B 1 (S n , m+1 + 5 n ,m-l)}, (C36) 

d^™ = A 2 sin (k±a)6 nm , (C37) 


which can be used similarly to generate the Floquet Green’s functions. 
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